#include <arm_neon.h>
#include <pthread.h>
#include <stdlib.h>
#include <iostream>
#include <sys/time.h>
#include <semaphore.h>
#include <unistd.h>
#include <omp.h>
using namespace std;
void OpenMP_Static_SIMD(float** a, int n)
{
#pragma omp parallel num_threads(8) shared(a,n)
    {
        for (int k = 0; k < n; k++)
        {
#pragma omp for simd schedule(static)
            for (int j = k + 1; j < n; j++)
                a[k][j] = a[k][j] / a[k][k];
#pragma omp barrier
            ;
#pragma omp critical
            a[k][k] = 1.0;
#pragma omp barrier
            ;
#pragma omp for simd collapse(2) schedule(static)
            for (int i = k + 1; i < n; i++)
                for (int j = k + 1; j < n; j++)
                    a[i][j] = a[i][j] - a[i][k] * a[k][j];
#pragma omp barrier
            ;
#pragma omp for simd schedule(static)
            for (int i = k + 1; i < n; i++)
                a[i][k] = 0;
        }
    }
}
void OpenMP_Static_Parallel(float** a, int n)
{
#pragma omp parallel for schedule(static) num_threads(8) shared(a,n)
    for (int k = 0; k < n; k++)
    {
        int j;
        for (j = k + 1; j % 4 != 0 && j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        for (; j <= n - 4; j += 4)
        {
            float32x4_t temp;
            float32x4_t t = vmovq_n_f32(a[k][k]);
            temp = vld1q_f32(a[k] + j);
            temp = vdivq_f32(temp, t);
            vst1q_f32(a[k] + j, temp);
        }
        for (; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            int j;
            for (j = k + 1; j % 4 != 0 && j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            for (; j <= n - 4; j += 4)
            {
                float32x4_t t1;
                float32x4_t t2;
                t1 = vld1q_f32(a[i] + j);
                t2 = vld1q_f32(a[k] + j);
                float32x4_t t3 = vmovq_n_f32(a[i][k]);
                t3 = vmulq_f32(t3, t2);
                t1 = vsubq_f32(t1, t3);
                vst1q_f32(a[i] + j, t1);
            }
            for (; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
        }
        for (int i = k + 1; i < n; i++)
            a[i][k] = 0;
    }
}
void OpenMP_Dynamic_SIMD(float** a, int n)
{
#pragma omp parallel for simd num_threads(8) shared(a,n)
    for (int k = 0; k < n; k++)
    {
        for (int j = k + 1; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            for (int j = k + 1; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            a[i][k] = 0;
        }
    }
}
void OpenMP_Dynamic_Parallel(float** a, int n)
{
#pragma omp parallel for schedule(dynamic) num_threads(8) shared(a,n)
    for (int k = 0; k < n; k++)
    {
        int j;
        for (j = k + 1; j % 4 != 0 && j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        for (; j <= n - 4; j += 4)
        {
            float32x4_t temp;
            float32x4_t t = vmovq_n_f32(a[k][k]);
            temp = vld1q_f32(a[k] + j);
            temp = vdivq_f32(temp, t);
            vst1q_f32(a[k] + j, temp);
        }
        for (; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            int j;
            for (j = k + 1; j % 4 != 0 && j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            for (; j <= n - 4; j += 4)
            {
                float32x4_t t1;
                float32x4_t t2;
                t1 = vld1q_f32(a[i] + j);
                t2 = vld1q_f32(a[k] + j);
                float32x4_t t3 = vmovq_n_f32(a[i][k]);
                t3 = vmulq_f32(t3, t2);
                t1 = vsubq_f32(t1, t3);
                vst1q_f32(a[i] + j, t1);
            }
            for (; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
        }
        for (int i = k + 1; i < n; i++)
            a[i][k] = 0;
    }
}
void OpenMP_Guided_SIMD(float** a, int n)
{
#pragma omp parallel num_threads(8) shared(a,n)
    {
        for (int k = 0; k < n; k++)
        {
#pragma omp for simd schedule(guided)
            for (int j = k + 1; j < n; j++)
                a[k][j] = a[k][j] / a[k][k];
#pragma omp barrier
            ;
#pragma omp critical
            a[k][k] = 1.0;
#pragma omp barrier
            ;
#pragma omp for simd collapse(2) schedule(guided)
            for (int i = k + 1; i < n; i++)
                for (int j = k + 1; j < n; j++)
                    a[i][j] = a[i][j] - a[i][k] * a[k][j];
#pragma omp barrier
            ;
#pragma omp for simd schedule(guided)
            for (int i = k + 1; i < n; i++)
                a[i][k] = 0;
        }
    }
}
void OpenMP_Guided_Parallel(float** a, int n)
{
#pragma omp parallel for schedule(guided) num_threads(8) shared(a,n)
    for (int k = 0; k < n; k++)
    {
        int j;
        for (j = k + 1; j % 4 != 0 && j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        for (; j <= n - 4; j += 4)
        {
            float32x4_t temp;
            float32x4_t t = vmovq_n_f32(a[k][k]);
            temp = vld1q_f32(a[k] + j);
            temp = vdivq_f32(temp, t);
            vst1q_f32(a[k] + j, temp);
        }
        for (; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            int j;
            for (j = k + 1; j % 4 != 0 && j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            for (; j <= n - 4; j += 4)
            {
                float32x4_t t1;
                float32x4_t t2;
                t1 = vld1q_f32(a[i] + j);
                t2 = vld1q_f32(a[k] + j);
                float32x4_t t3 = vmovq_n_f32(a[i][k]);
                t3 = vmulq_f32(t3, t2);
                t1 = vsubq_f32(t1, t3);
                vst1q_f32(a[i] + j, t1);
            }
            for (; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
        }
        for (int i = k + 1; i < n; i++)
            a[i][k] = 0;
    }
}
void m_rest(float** a, int n)
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < i; j++)
            a[i][j] = 0;
        a[i][i] = 1.0;
        for (int j = i + 1; j < i; j++)
        {
            float x = rand();
            while (x == 0) x = rand();
            a[i][j] = x;
        }
    }
    for (int k = 0; k < n; k++)
        for (int i = k + 1; i < n; i++)
            for (int j = 0; j < n; j++)
                a[i][j] += a[k][j];
}
int main()
{
    srand((unsigned)time(NULL));
    int n = 10;
    while (n < 1600)
    {
        float** a = (float**)aligned_alloc(32, n * sizeof(float*));
        for (int i = 0; i < n; i++)
            a[i] = (float*)aligned_alloc(32, n * sizeof(float));
        m_rest(a, n);

        cout << n << " ";
        struct timeval tstart, tend;
        double timeUsed = 0.0;
        float** t = (float**)aligned_alloc(32, n * sizeof(float*));
        for (int i = 0; i < n; i++)
            t[i] = (float*)aligned_alloc(32, n * sizeof(float));


        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        OpenMP_Static_SIMD(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        cout << timeUsed << " ";

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        OpenMP_Static_Parallel(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        cout << timeUsed << " ";

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        OpenMP_Dynamic_SIMD(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        cout << timeUsed << " ";

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        OpenMP_Dynamic_Parallel(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        cout << timeUsed << " ";

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        OpenMP_Guided_SIMD(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        cout << timeUsed << " ";

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        OpenMP_Guided_Parallel(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        cout << timeUsed << " ";

        cout << endl;
        for (int i = 0; i < n; i++)
        {
            delete[] t[i];
            delete[] a[i];
        }
        delete[] t;
        delete[] a;
        if (n >= 800) n += 200;
        if (n >= 400 && n < 800) n += 100;
        if (n >= 200 && n < 400) n += 50;
        if (n < 200) n += 10;
    }
    return 0;
}
